/*------------------------------------------------------------------------
 *  Copyright (C) 2010  Luis M. de la Cruz
 *
 *  This file is part of TUNA
 *
 *  TUNA is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  TUNA is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ------------------------------------------------------------------------*/

namespace Tuna {

  
/*-----------------------------------------------------------------------
 * DESCRIPTION
 *     General algorithm for to solve tridiagonal systems like this:
 *
 *     |a1   b1    0    0    0  ...    0 |   | x1 |   | d1 |
 *     |c2   a2   b2    0    0  ...    0 |   | x2 |   | d2 |
 *     | 0   c3   a3   b3    0  ...    0 |   | x3 |   | d3 |
 *     | 0    0   c4   a4   b4  ...    0 | * | x4 | = | d4 |
 *     | .    .    .    .    .  ...    0 |   |  . |   |  . |
 *     | .    .    .    .    .  ...    0 |   |  . |   |  . |
 *     | 0    0    0    0  cN-1 aN-1 bN-1|   |xN-1|   |dN-1|
 *     | 0    0    0    0    0   cN   aN |   | xN |   | dN |
 *  
 *     I'm using the Thomas Algorithm:
 *     
 *     P1 = b1 / a1; Q1 = d1 / a1
 *
 *     for k = 2 to N do              // Forward step
 *        Pk = bk / ( ak - ck * Pk-1 )  
 *        Qk = ( dk - ck * Qk-1 ) / ( ak - ck * Pk-1 )
 *
 *     xN = QN
 *
 *     for k = N-1 to 1               // Backward substitution
 *        xk = Qk - Pk * xk+1 
 *    
 *     x contains the final solution.
 *
 *-----------------------------------------------------------------------*/

template<typename Matrix, typename Vector>
bool
Solver::thomas(const Matrix &A, Vector &x, const Vector &d)
{
  int N = x.size(); 
  
  Vector P(N), Q(N);
  const Vector &b = A.Storage.b;
  const Vector &a = A.Storage.a;
  const Vector &c = A.Storage.c;
  
  std::cout << a << std::endl;
  std::cout << b << std::endl;
  std::cout << c << std::endl;

  // In this implementation arrays begin in 0
  P(0) = b(0) / a(0);
  Q(0) = d(0) / a(0);
  
  // Forward step 
  double divisor;
  
  for(int k = 1; k < N; ++k)
    {
      divisor = a(k) - c(k) * P(k-1);
      P(k) = b(k) / divisor;
      Q(k) = ( d(k) - c(k) * Q(k-1) ) / divisor;
    }

  x(N-1) = Q(N-1);
  
  // Backward substitution 
  for(int k = N-2; k >= 0; --k)
    x(k) = Q(k) - P(k) * x(k+1);
  
  return 0;
}

} // Tuna namespace



